package scu.maqiang.numeric;

import java.util.Arrays;
import java.util.Random;
import java.util.Scanner;

import scu.maqiang.numeric.Pair;
import java.util.stream.DoubleStream;
import java.util.stream.IntStream;

import static java.lang.Math.*;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileWriter;
import java.io.IOException;

/**
 * 矩阵与向量操作静态类
 * 
 * @author 马强
 *
 */
public class MVO {
	
	/**
	 * 生成一个的随机向量
	 * @param m 向量长度
	 * @return 长度为m的随机向量
	 */
	public static double[] rand(int m) {
		double[] result = new double[m];
		Random r = new Random();
		Arrays.setAll(result, i -> r.nextDouble());
		return result;
	}

	/**
	 * 生成一个0矩阵
	 * @param m
	 * @param n
	 * @return
	 */
	public static double[][] zero(int m, int n) {
		double[][] result = new double[m][];
		for(int i = 0; i < m; ++i) {
			result[i] = new double[n];
			Arrays.fill(result[i], 0.0);
		}
		return result;
	}


	/**
	 * 生成一个随机矩阵
	 * @param m 矩阵行数
	 * @param n 矩阵列数
	 * @return 行数为m，列数为n的随机矩阵
	 */
	public static double[][] rand(int m, int n) {
		double[][] result = new double[m][];
		Random r = new Random();
		for(int i = 0; i < m; i++) {
			result[i] = new double[n];
			Arrays.setAll(result[i], j -> r.nextDouble());
		}
		return result;
	}
	
	/**
	 * 生成一个单位矩阵
	 * @param n 矩阵阶数
	 * @return n阶单位矩阵
	 */
	public static double[][] eye(int n) {
		double[][] result = new double[n][n];
		for(int i = 0; i < n; i++) {
			result[i][i] = 1.0;
		}
		return result;
	}

	/**
	 * 生成一个向量, 设置向量元素为常值
	 * @param n 向量长度
	 * @param value 向量元素的值
	 * @return 长度为n, 各元素均为value的向量
	 */
	public static double[] vec(int n, double value) {
		double[] result = new double[n];
		Arrays.fill(result, value);
		return result;
	}
	/**
	 * 计算一串实数的最大值, 如果没有实数，则返回NAN
	 * @param x 实数可变参数
	 * @return 一串实数中的最大值
	 */
	public static double max(double... x) {
		return DoubleStream.of(x).max().orElse(Double.NaN);
	}

	/**
	 * 计算一串实数的最小值, 如果没有实数，则返回NANֵ
	 * @param x 实数可变参数
	 * @return 一串实数中的最小值ֵ
	 */
	public static double min(double... x) {
		return DoubleStream.of(x).min().orElse(Double.NaN);
	}

	/**
	 * 得到二维数组的某一行所有元素
	 * 
	 * @param a 二维矩阵数组
	 * @param i 矩阵a的第i行, 要求大于等0， 小于矩阵的行数
	 * @return 矩阵a的第i行元素
	 */
	public static double[] row(double[][] a, int i) {
		ParamCheck.checkBounds(i, 0, a.length);
		return a[i].clone();
	}

	/**
	 *  得到二维数组的某一列所有元素
	 * 
	 * @param a 二维矩阵数组
	 * @param i 矩阵a的第i列, 要求大于等0， 小于矩阵的列数
	 * @return 矩阵a的第i列元素
	 */
	public static double[] col(double[][] a, int i) {
		ParamCheck.checkBounds(i, 0, a[0].length);
		double[] newCol = new double[a.length];
		for (int j = 0; j < a.length; j++) {
			newCol[j] = a[j][i];
		}
		return newCol;
	}

	public static double[] constant(int size, double val) {
		double[] result = new double[size];
		Arrays.fill(result, val);
		return result;
	}

	/**
	 * 计算一串整数的最大值, 如果没有实数，则返回MAX_VALUE
	 * 
	 * @param x 整数可变参数
	 * @return 一串整数中的最大值ֵ
	 */
	public static int max(int... x) {
		return IntStream.of(x).max().orElse(Integer.MAX_VALUE);
	}

	/**
	 * 得到一列实数中绝对值最大值
	 * 
	 * @param x
	 *            可变实数数组
	 * @return 实数中绝对值最大值，如果未传入实数，返回NAN
	 */
	public static double maxA(double... x) {
		return DoubleStream.of(x).map(Math::abs).max().orElse(Double.NaN);
	}


	/**
	 * 得到实数数组中的最大值以及相应的位置
	 * 
	 * @param x
	 *            双精度数组
	 * @param pair
	 *            整数与实数组成的数对，分别得到实数数组中最大值及其所在位置
	 */
	public static void max(double[] x, Pair<Integer, Double> pair) {
		double maxval = x[0];
		int maxidx = 0;
		for (int i = 1, n = x.length; i < n; i++) {
			if (maxval < x[i]) {
				maxval = x[i];
				maxidx = i;
			}
		}
		pair.setFirst(maxidx);
		pair.setSecond(maxval);
	}

	/**
	 * 求整数数组的最大值以及最大值所在索引
	 * @param x 整数数组
	 * @param pair 数对对象, 第一个值为数组最大值, 第二个值为最大值所在数组的索引
	 */
	public static void max(int[] x, Pair<Integer, Integer> pair) {
		int maxval = x[0];
		int maxidx = 0;
		for (int i = 1, n = x.length; i < n; i++) {
			if (maxval < x[i]) {
				maxval = x[i];
				maxidx = i;
			}
		}
		pair.setFirst(maxidx);
		pair.setSecond(maxval);
	}

	public static Pair<Integer, Double> maxAbsColValDown(double[][] A, int i) {
		int n = A.length;
		ParamCheck.checkBounds(i, 0, n);
		double maxval = A[i][i];
		int maxIdx = i;
		for (int idx = i + 1; idx < n; idx++) {
			if (Math.abs(A[idx][i]) > Math.abs(maxval)) {
				maxval = A[idx][i];
				maxIdx = idx;
			}
		}
		return new Pair<>(maxIdx, maxval);
	}

	/**
	 * 求实数数组的平均值
	 * 
	 * @param x
	 *            可变实数数组
	 * @return 数组元素的平均值
	 */
	public static double average(double... x) {
		return DoubleStream.of(x).average().orElse(Double.NaN);
	}

	/**
	 * 求整数数组的平均值
	 * 
	 * @param x
	 *            可变整数数组
	 * @return 数组元素的平均值
	 */
	public static double average(int... x) {
		return IntStream.of(x).average().orElse(Double.NaN);
	}

	/**
	 * 求数组元素的最小值以及相应的下标
	 * 
	 * @param x
	 *            实数数组
	 * @param pair
	 *            数对，第一个值为最小值下标，第二个值为最小值
	 */
	public static void min(double[] x, Pair<Integer, Double> pair) {
		double minval = x[0];
		int minidx = 0;
		for (int i = 1, n = x.length; i < n; i++) {
			if (minval > x[i]) {
				minval = x[i];
				minidx = i;
			}
		}
		pair.setFirst(minidx);
		pair.setSecond(minval);
	}

	public static void min(int[] x, Pair<Integer, Integer> pair) {
		int minval = x[0];
		int minidx = 0;
		for (int i = 1, n = x.length; i < n; i++) {
			if (minval > x[i]) {
				minval = x[i];
				minidx = i;
			}
		}
		pair.setFirst(minidx);
		pair.setSecond(minval);
	}

	public static int min(int... x) {
		return IntStream.of(x).min().orElse(Integer.MAX_VALUE);
	}

	public static boolean isZero(double[] x) {
		for(int i = 0, n = x.length; i < n; ++i) {
			if (x[i] != 0) {
				return false;
			}
		}
		return true;
	}

	public static boolean isZero(double[][] x) {
		for(int i = 0, n = x.length; i < n; ++i) {
			for(int j = 0, nj = x[i].length; j < nj; ++j) {
				if (x[i][j] != 0) {
					return false;
				}
			}
		}
		return true;
	}
	/**
	 * 计算实数数组和
	 * 
	 * @param x
	 *            可变实数数组
	 * @return 实数数组各元素之和
	 */
	public static double sum(double... x) {
		return DoubleStream.of(x).sum();
	}

	public static Complex sum(Complex... x) {
		Complex result = new Complex(0.0, 0.0);
		for(int i = 0, n = x.length; i < n; i++) {
			result.real = x[i].real;
			result.image = x[i].image;
		}
		return result;
	}

	public static int sum(int... x) {
		return IntStream.of(x).sum();
	}

	public static double[] sumRow(double[][] x) {
		double[] result = new double[x.length];
		Arrays.setAll(result, i -> DoubleStream.of(x[i]).sum());
		return result;
	}

	public static double[] sumCol(double[][] x) {
		int m = x.length, n = x[0].length;
		double[] result = new double[n];
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < m; j++) {
				result[i] += x[j][i];
			}
		}
		return result;
	}

	public static double[] getSlice(double[] u, int[] idx) {
		int n = idx.length;
		double[] result = new double[n];
		for (int i = 0; i < n; i++) {
			result[i] = u[idx[i]];
		}
		return result;
	}

	/**
	 * 
	 * @param u
	 * @param idx
	 * @param result
	 */
	public static void getSlice(double[] u, int[] idx, double[] result) {
		int n = idx.length;
		ParamCheck.checkEqual(n, result.length);
		for (int i = 0; i < n; i++) {
			result[i] = u[idx[i]];
		}
	}

	/**
	 * 计算实数数组的L2范数
	 * 
	 * @param x
	 *            可变实数数组
	 * @return L2范数，各元素平方和再开方取正数
	 */
	public static double L2Norm(double... x) {
		// return sqrt(DoubleStream.of(x).map(ele -> ele * ele).sum());
		double temp = 0.0;
		for (int i = 0, n = x.length; i < n; i++) {
			temp += x[i] * x[i];
		}
		return sqrt(temp);
	}
	
	public static double L2Norm2(double... x) {
		double result = 0.0;
		for (int i = 0, n = x.length; i < n; i++) {
			result += x[i] * x[i];
		}
		return result;		
	}

	/**
	 * 计算数组x1加上一个数alpha乘以x2之和，即计算x1 + alpha * x2
	 * 
	 * @param x1
	 *            实数数组
	 * @param alpha
	 *            实数
	 * @param x2
	 *            实数数组
	 * @return 实数数组，结果为：x1 + alpha * x2， 当alpha = 1时，可用于计算两数组之和，alpha=-1时，可用于计算两数组之差
	 */
	public static double[] add(double[] x1, double alpha, double[] x2) {
		int n = x1.length;
		double[] result = new double[n];
		Arrays.setAll(result, i -> x1[i] + alpha * x2[i]);
		return result;
	}

	public static void add(double[] x1, double alpha, double[] x2, double[] result) {
		int n = x1.length;
		ParamCheck.checkEqual(n, result.length);
		Arrays.setAll(result, i -> x1[i] + alpha * x2[i]);
	}

	/**
	 * 求两数组内积
	 * 
	 * @param x1
	 *            数组1
	 * @param x2
	 *            数组2
	 * @return 两数组内积，两数组各分量之积之和
	 */
	public static double dot_product(double[] x1, double[] x2) {
		ParamCheck.checkEqual(x1.length, x2.length);
		double result = 0.0;
		for (int i = 0, n = x1.length; i < n; i++) {
			result += x1[i] * x2[i];
		}
		return result;
	}

	/**
	 * 
	 * @param x1
	 * @param idx
	 * @param x2
	 * @return
	 */
	public static double dot_product(double[] x1, int[] idx, double[] x2) {
		int n = idx.length;
		ParamCheck.checkEqual(n, x2.length);
		double result = 0.0;
		for (int i = 0; i < n; i++) {
			result += x1[idx[i]] * x2[i];
		}
		return result;
	}

	public static double[] cross_product(double[] x, double[] y) {
		double[] result = new double[3];
		result[0] = x[1] * y[2] - x[2] * y[1];
		result[1] = x[2] * y[0] - x[0] * y[2];
		result[2] = x[0] * y[1] - x[1] * y[0];
		return result;
	}

	/**
	 * 
	 * @param x1
	 * @param alpha
	 * @param x2
	 * @return
	 */
	public static double[][] add(double[][] x1, double alpha, double[][] x2) {
		int m = x1.length;
		int n = x1[0].length;
		double[][] result = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				result[i][j] = x1[i][j] + alpha * x2[i][j];
			}
		}
		return result;
	}

	/**
	 * 
	 * @param z0
	 * @param d
	 * @return
	 */
	public static double[] add(double[] z0, double d) {
		int m = z0.length;
		double[] result = new double[m];
		for (int i = 0; i < m; i++) {
			result[i] = z0[i] + d;
		}
		return result;
	}

	public static void mulSelf(double[][] z0, double d) {
		int m = z0.length;
		int n = z0[0].length;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				z0[i][j] *= d;
			}
		}
	}

	/**
	 * 
	 * @param z0
	 * @param d
	 * @return
	 */
	public static double[][] mul(double[][] z0, double d) {
		int m = z0.length;
		int n = z0[0].length;
		double[][] result = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				result[i][j] = z0[i][j] * d;
			}
		}
		return result;
	}

	public static double[] mul(double[] r1, double[] r2) {
		int n = r1.length;
		double[] result = new double[n];
		ParamCheck.checkEqual(n, r2.length);
		for (int i = 0; i < n; i++) {
			result[i] = r1[i] * r2[i];
		}
		return result;
	}

	public static void mul(double[] r1, double[] r2, double[] result) {
		int n = r1.length;
		ParamCheck.checkEqual(n, r2.length);
		ParamCheck.checkEqual(n, result.length);
		for (int i = 0; i < n; i++) {
			result[i] = r1[i] * r2[i];
		}
	}

	public static void mul(double[][] z0, double d, double[][] result) {
		int m = z0.length;
		int n = z0[0].length;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				result[i][j] = z0[i][j] * d;
			}
		}
	}

	public static double[][] matmul(double[][] A, double[][] B) {
		int m = A.length;
		int k = A[0].length;
		ParamCheck.checkEqual(k, B.length);
		int n = B[0].length;
		double[][] C = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				double temp = 0.0;
				for (int kk = 0; kk < k; kk++) {
					temp += A[i][kk] * B[kk][j];
				}
				C[i][j] = temp;
			}
		}
		return C;
	}

	public static void matmul(double[][] A, double[][] B, double[][] C) {
		int m = A.length;
		int k = A[0].length;
		ParamCheck.checkEqual(k, B.length);
		int n = B[0].length;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				double temp = 0.0;
				for (int kk = 0; kk < k; kk++) {
					temp += A[i][kk] * B[kk][j];
				}
				C[i][j] = temp;
			}
		}
	}

	public static double[] matmul(double[][] A, double[] B) {
		int m = A.length;
		int n = A[0].length;
		ParamCheck.checkEqual(n, B.length);
		double[] C = new double[m];
		for (int i = 0; i < m; i++) {
			double temp = 0.0;
			for (int j = 0; j < n; j++) {
				temp += A[i][j] * B[j];
			}
			C[i] = temp;
		}
		return C;
	}

	public static void matmul(double[][] A, double[] B, double[] C) {
		int m = A.length;
		int n = A[0].length;
		ParamCheck.checkEqual(n, B.length);
		ParamCheck.checkEqual(m, C.length);
		for (int i = 0; i < m; i++) {
			double temp = 0.0;
			for (int j = 0; j < n; j++) {
				temp += A[i][j] * B[j];
			}
			C[i] = temp;
		}
	}

	public static double[] matmulT(double[][] A, double[] B) {
		int m = A.length;
		int n = A[0].length;
		ParamCheck.checkEqual(m, B.length);
		double[] result = new double[n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				result[j] += A[i][j] * B[i];
			}
		}
		return result;
	}

	public static void matmulT(double[][] A, double[] B, double[] C) {
		int m = A.length;
		int n = A[0].length;
		ParamCheck.checkEqual(m, B.length);
		ParamCheck.checkEqual(n, C.length);
		Arrays.fill(C, 0.0);
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				C[j] += A[i][j] * B[i];
			}
		}
	}

	public static Complex[] matmul(Complex[][] mat, Complex[] vec) {
		int m = mat.length;
		int n = mat[0].length;
		ParamCheck.checkEqual(m, vec.length);
		Complex[] result = new Complex[m];
		for (int i = 0; i < m; i++) {
			double re = 0.0;
			double im = 0.0;
			for (int j = 0; j < n; j++) {
				re += mat[i][j].real * vec[j].real - mat[i][j].image * vec[j].image;
				im += mat[i][j].real * vec[j].image + mat[i][j].image * vec[j].real;
			}
			vec[i].real = re;
			vec[i].image = im;
		}
		return result;
	}

	public static double[][] transpose(double[][] A) {
		int m = A.length;
		int n = A[0].length;
		double[][] At = new double[n][m];
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < m; j++) {
				At[i][j] = A[j][i];
			}
		}
		return At;
	}

	public static void transpose(double[][] A, double[][] At) {
		int m = A.length;
		int n = A[0].length;
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < m; j++) {
				At[i][j] = A[j][i];
			}
		}
	}

	public static void setValue(double[] a, double val) {
		Arrays.fill(a, val);
	}

	public static void fill(double[][] a, double value) {
		for (double[] doubles : a) {
			Arrays.fill(doubles, value);
		}
	}

	public static double[] linspace(double a, double b, int n) {

		double[] v = new double[n];
		double dx = (b - a) / (n - 1);
		Arrays.setAll(v, i -> a + i * dx);
		return v;
	}

	public static double distance(double[] x0, double[] x1) {
		ParamCheck.checkEqual(x0.length, x1.length);
		double result = 0.0;
		double temp;
		for (int i = 0, n = x0.length; i < n; i++) {
			temp = x0[i] - x1[i];
			result += temp * temp;
		}
		return sqrt(result);
	}

	public static double yAx(double[] y, double[][] A, double[] x) {
		int m = y.length;
		int n = x.length;
		ParamCheck.checkEqual(m, A.length);
		ParamCheck.checkEqual(A[0].length, n);
		double result = 0.0;
		double temp;
		for (int i = 0; i < m; i++) {
			temp = 0.0;
			for (int j = 0; j < n; j++) {
				temp += A[i][j] * x[j];
			}
			result += y[i] * temp;
		}
		return result;
	}

	public static double determinant(double[][] A) {
		if (A.length == 4) {
			double[][] M11 = { { A[1][1], A[1][2], A[1][3] }, { A[2][1], A[2][2], A[2][3] },
					{ A[3][1], A[3][2], A[3][3] } };
			double[][] M12 = { { A[1][0], A[1][2], A[1][3] }, { A[2][0], A[2][2], A[2][3] },
					{ A[3][0], A[3][2], A[3][3] } };
			double[][] M13 = { { A[1][0], A[1][1], A[1][3] }, { A[2][0], A[2][1], A[2][3] },
					{ A[3][0], A[3][1], A[3][3] } };
			double[][] M14 = { { A[1][0], A[1][1], A[1][2] }, { A[2][0], A[2][1], A[2][2] },
					{ A[3][0], A[3][1], A[3][2] } };
			return A[0][0] * determinant(M11) - A[0][1] * determinant(M12) + A[0][2] * determinant(M13)
					- A[0][3] * determinant(M14);
		}
		if (A.length == 3) {
			return A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - A[0][1] * (A[1][0] * A[2][2] - A[2][0] * A[1][2])
					+ A[0][2] * (A[1][0] * A[2][1] - A[2][0] * A[1][1]);
		}
		if (A.length == 2) {
			return A[0][0] * A[1][1] - A[0][1] * A[1][0];
		}
		if (A.length == 1) {
			return A[0][0];
		}
		throw new IllegalArgumentException("The size of matrix should not be greater than 4!");
	}


	public static double[] getThickness(double[][] boxMatrix) {
		double[] thickness = new double[3];
		double volume = Math.abs(MVO.determinant(boxMatrix));
		double[] a = {boxMatrix[0][0], boxMatrix[1][0], boxMatrix[2][0]};
		double[] b = {boxMatrix[0][1], boxMatrix[1][1], boxMatrix[2][1]};
		double[] c = {boxMatrix[0][2], boxMatrix[1][2], boxMatrix[2][2]};
		thickness[0] = volume / MVO.L2Norm(MVO.cross_product(b, c));
		thickness[1] = volume / MVO.L2Norm(MVO.cross_product(c, a));
		thickness[2] = volume / MVO.L2Norm(MVO.cross_product(a, b));
		return thickness;
	}

	public static double[] getThickness2D(double[][] boxMatrix) {
		double[] thickness = new double[2];
		double area = abs(boxMatrix[0][0] * boxMatrix[1][1] - boxMatrix[0][1] * boxMatrix[1][0]);
		double lenX = sqrt(boxMatrix[0][0] * boxMatrix[0][0] + boxMatrix[1][0] * boxMatrix[1][0]);
		double lenY = sqrt(boxMatrix[0][1] * boxMatrix[0][1] + boxMatrix[1][1] * boxMatrix[1][1]);
		thickness[0] = area / lenY;
		thickness[1] = area / lenX;
		return thickness;
	}


	public static String toString(double... v) {
		StringBuilder sb = new StringBuilder();
		sb.append(v.length + "\n");
		DoubleStream.of(v).forEach(ele -> sb.append(String.format("%30.20g\n", ele)));
		return sb.append("\n").toString();
	}

	/**
	 * 
	 * @param a
	 * @return
	 */
	public static String toString(double[][] a) {
		StringBuffer sb = new StringBuffer("");
		for (double[] row : a) {
			for (double ele : row) {
				sb.append(String.format("%20.10e", ele));
			}
			sb.append("\n");
		}
		return sb.toString();
	}

	public static String toString(int[][] a) {
		System.out.println(a.length + "\t  " + a[0].length);
		StringBuffer sb = new StringBuffer("");
		for (int[] row : a) {
			for (int ele : row) {
				sb.append(String.format("%10d", ele));
			}
			sb.append("\n");
		}
		return sb.toString();
	}

	/**
	 * 以二维数组为矩阵，提出取矩阵相应子矩阵
	 * @param A 原矩阵
	 * @param row0 行开始指标，下标从0记数
	 * @param row1 行结束指标，下标从0记数
	 * @param col0 列开始指标，下标从0记数
	 * @param col1 列结束指标，下标从0记数
	 * @return
	 */
	public static double[][] subMatrix(double[][] A, int row0, int row1, int col0, int col1) {
		int m = row1 - row0 + 1;
		int n = col1 - col0 + 1;
		double[][] result = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				result[i][j] = A[row0 + i][col0 + j];
			}
		}
		return result;
	}

	public static double[] subVector(double[] v, int start, int end) {
		int n = end - start + 1;
		double[] result = new double[n];
		for (int i = 0; i < n; i++) {
			result[i] = v[i + start];
		}
		return result;
	}

	/**
	 * 
	 * @param A
	 * @param B
	 */
	public static void matmul_toRight(double[][] A, double[][] B) {
		int m = A.length;
		int n = B[0].length;
		ParamCheck.checkEqual(m, A[0].length);
		ParamCheck.checkEqual(m, B.length);
		double[][] temp = new double[m][n];

		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				double tt = 0.0;
				for (int k = 0; k < m; k++) {
					tt += A[i][k] * B[k][j];
				}
				temp[i][j] = tt;
			}
		}
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				B[i][j] = temp[i][j];
			}
		}
	}

	/**
	 * 
	 * @param A
	 * @param B
	 */
	public static void matmul_toLeft(double[][] A, double[][] B) {
		int m = A.length;
		int n = B[0].length;
		ParamCheck.checkEqual(n, A[0].length);
		ParamCheck.checkEqual(n, B.length);
		double[][] temp = new double[m][n];
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				double tt = 0.0;
				for (int k = 0; k < n; k++) {
					tt += A[i][k] * B[k][j];
				}
				temp[i][j] = tt;
			}
		}
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				A[i][j] = temp[i][j];
			}
		}
	}

	/**
	 * 使用Householder变换对矩阵A作QR分解
	 * @param A mxn型矩阵，要求m>n
	 * @param Q mxm型正交矩阵
	 * @param R mxn型矩阵
	 */
	public static void QR_HH(double[][] A, double[][] Q, double[][] R) {
		int m = A.length;
		int n = A[0].length;
		ParamCheck.checkEqual(m, Q.length);
		ParamCheck.checkEqual(m, Q[0].length);
		ParamCheck.checkEqual(m, R.length);
		ParamCheck.checkEqual(n, R[0].length);
		fill(Q, 0.0);
		for (int i = 0; i < m; i++) {
			Q[i][i] = 1.0;
		}
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				R[i][j] = A[i][j];
			}
		}
		double[][] H = new double[m][m];
		double[] w = new double[m];
		for (int i = 0; i < n; i++) {
			fill(H, 0.0);
			setValue(w, 0.0);
			for (int j = 0; j < i; j++) {
				H[j][j] = 1.0;
			}
			double[] coli = MVO.col(R, i);
			double temp = 0.0;
			for (int k = i; k < m; k++) {
				temp += coli[k] * coli[k];
			}
			w[i] = Math.sqrt(temp);
			double[] v = MVO.add(w, -1.0, coli);
			double L2V = 0.0;
			for (int k = i; k < m; k++) {
				L2V += v[k] * v[k];
			}
			for (int ii = i; ii < m; ii++) {
				for (int jj = i; jj < m; jj++) {
					double diag = ii == jj ? 1.0 : 0.0;
					H[ii][jj] = diag - 2 * v[ii] * v[jj] / L2V;
				}
			}
			matmul_toRight(H, R);
			matmul_toLeft(Q, H);
		}
	}

	public static double[][] UpperTriangular(int n) {
		double[][] ut = new double[n][];
		for (int i = 0; i < n; i++) {
			ut[i] = new double[n - i];
		}
		return ut;
	}

	public static double[][] LowerTriangular(int n) {
		double[][] lt = new double[n][];
		for (int i = 0; i < n; i++) {
			lt[i] = new double[i + 1];
		}
		return lt;
	}

	public static double[][] upperHesenberger(int m, int n) {
		double[][] uh = new double[m][];
		uh[0] = new double[n];
		for (int i = 1; i < m; i++) {
			uh[i] = new double[n - i + 1];
		}
		return uh;
	}

	public static int[] readIntArray(String fileName) {
		int[] result = null;
		try (Scanner sc = new Scanner(new File(fileName))) {
			int len = sc.nextInt();
			result = new int[len];
			for (int i = 0; i < len; i++) {
				result[i] = sc.nextInt();
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
		return result;
	}

	public static double[] readArray(String fileName) {
		double[] result = null;
		try (Scanner sc = new Scanner(new File(fileName))) {
			int len = sc.nextInt();
			result = new double[len];
			for (int i = 0; i < len; i++) {
				result[i] = sc.nextDouble();
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
		return result;
	}

	public static double[][] readArray2(String fileName) {
		double[][] result = null;
		try (Scanner sc = new Scanner(new File(fileName))) {
			int row = sc.nextInt();
			int col = sc.nextInt();
			result = new double[row][col];
			for (int i = 0; i < row; i++) {
				for(int j = 0; j < col; j++) {
					result[i][j] = sc.nextDouble();
				}
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
		return result;
	}

	public static void saveArray(String fileName, double[] vec) {
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			int n = vec.length;
			bw.write(n + "\n");
			for (int i = 0; i < n; i++) {
				bw.write(vec[i] + "\n");
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	public static void saveMatrix(String fileName, int[][] mat) {
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			int m = mat.length;
			int n = mat[0].length;
			bw.write(m + " \t" + n + "\n");
			for (int i = 0; i < m; i++) {
				for(int j = 0; j < n; j++) {
					bw.write(mat[i][j] + " \t");
				}
				bw.write("\n");
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	public static void saveMatrix(String fileName, double[][] mat) {
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			int m = mat.length;
			int n = mat[0].length;
			bw.write(m + " \t" + n + "\n");
			for (int i = 0; i < m; i++) {
				for(int j = 0; j < n; j++) {
					bw.write(mat[i][j] + " \t");
				}
				bw.write("\n");
			}
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

	public static double polyval(double[] coef, double[] b, double x) {
		int n = coef.length - 1;
		double result = coef[n];
		for (int i = n - 1; i >= 0; i--) {
			result = result * (x - b[i]) + coef[i];
		}
		return result;
	}

	public static double polyval(double[] coef, double x) {
		int n = coef.length - 1;
		double result = coef[n];
		for (int i = n - 1; i >= 0; i--) {
			result = result * x + coef[i];
		}
		return result;
	}

	public static double[] ployval(double[] coef, double[] b, double[] t) {
		int n = coef.length - 1;
		double result[] = new double[t.length];
		for(int k = 0; k < t.length; k++) {
			result[k] = coef[n];
			for (int i = n - 1; i >= 0; i--) {
				result[k] = result[k] * (t[k] - b[i]) + coef[i];
			}
		}
		return result;
	}

	public static double[] ployval(double[] coef, double[] t) {
		int n = coef.length - 1;
		double result[] = new double[t.length];
		for(int k = 0; k < t.length; k++) {
			result[k] = coef[n];
			for (int i = n - 1; i >= 0; i--) {
				result[k] = result[k] * t[k] + coef[i];
			}
		}
		return result;
	}


	public static void toFile(double[] x, double[] y, String fileName) {
		int n = x.length;
		ParamCheck.checkEqual(n, y.length);
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			for (int i = 0; i < n; i++) {
				bw.write(x[i] + "\t" + y[i] + "\n");
			}
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	}

	public static void toFile(double[] x, double[][] y, String fileName) {
		int n = x.length;
		int m = y.length;
		for (int j = 0; j < m; j++) {
			ParamCheck.checkEqual(n, y[j].length);
		}
		try (BufferedWriter bw = new BufferedWriter(new FileWriter(fileName))) {
			for (int i = 0; i < n; i++) {
				bw.write(x[i] + "\t");
				for (int j = 0; j < m; j++) {
					bw.write(y[j][i] + "\t");
				}
				bw.write("\n");
			}
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	}
	
	public static boolean match(int[] inci1, int[] inci2) {
		int n = inci1.length;
		ParamCheck.checkEqual(n, inci2.length);
		boolean[] flag = new boolean[n];
		int count = 0;
		for(int i = 0; i < n; i++) {
			for(int j = 0; j < n; j++) {
				if(!flag[j] && inci1[i] == inci2[j]) {
					flag[j] = true;
					count++;
					break;
				}
			}
		}
		return n == count;
	}
	
	public static void main(String[] args) {
		int[] idx1 = {1, 2, 5};
		int[] idx2 = {2, 1, 4};
		System.out.println(MVO.match(idx1, idx2));
		idx2[2] = 5;
		System.out.println(MVO.match(idx1, idx2));
		
		int[] idx3 = {3, 5, 7, 8};
		int[] idx4 = {5, 3, 8, 7};
		System.out.println(MVO.match(idx3, idx4));
		int[] idx5 = {5, 3, 8, 6};
		System.out.println(MVO.match(idx3, idx5));
	}

	public static void fill(Complex[] v, double value) {
		int n = v.length;
		for(int i = 0; i < n; i++) {
			v[i].real = value;
			v[i].image = 0.0;
		}
	}
	
	public static double[] removeIdx(double[] vec, int idx) {
		int n = vec.length;
		ParamCheck.isPositive(n);
		ParamCheck.checkBounds(idx, 0, n);
		int nn = n - 1;
		double[] result = new double[nn];
		int count = 0;
		for(int i = 0; i < n; i++) {
			if(i != idx) {
				result[count++] = vec[i];
			}
		}
		return result;
	}

	public static double[][] inverseMat2(double[][] mat) {
		double det = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
		if(Math.abs(det) > 1.0e-10) {
			return new double[][]{{mat[1][1]/det, -mat[0][1]/det}, {-mat[1][0]/det, mat[0][0]/det}};
		} else {
			System.out.println("Matrix cannot be inversed!");
			return null;
		}
	}

	public static double[][] inverseMat3(double[][] box) {
		double[][] invBox = MVO.zero(3, 3);
		invBox[0][0] = box[1][1] * box[2][2] - box[1][2] * box[2][1];
		invBox[0][1] = box[0][2] * box[2][1] - box[0][1] * box[2][2];
		invBox[0][2] = box[0][1] * box[1][2] - box[0][2] * box[1][1];
		invBox[1][0] = box[1][2] * box[2][0] - box[1][0] * box[2][2];
		invBox[1][1] = box[0][0] * box[2][2] - box[0][2] * box[2][0];
		invBox[1][2] = box[0][2] * box[1][0] - box[0][0] * box[1][2];
		invBox[2][0] = box[1][0] * box[2][1] - box[1][1] * box[2][0];
		invBox[2][1] = box[1][0] * box[2][0] - box[0][0] * box[2][1];
		invBox[2][2] = box[0][0] * box[1][1] - box[0][1] * box[1][0];
		double det = MVO.determinant(box);
		for (int i = 0; i < 3; ++i) {
			for (int j = 0; j < 3; ++j) {
				invBox[i][j] /= det;
			}
		}
		return invBox;
	}

	/**
	 * 矩阵A的不完全Cholesky分解
	 * @param A
	 * @return
	 */
	public static double[][] iCholesky(double[][] A) {
		int m = A.length;
		int n = A[0].length;
		assert (m == n);
		double[][] R = new double[n][n];
		for(int i = 0; i < n; i++) {
			double temp = 0.0;
			for(int j = 0; j <= i - 1; j++) {
				temp += R[j][i] * R[j][i];
			}
			double temp2 = A[i][i] - temp;
			if (temp2 <= 0) {
				throw new IllegalArgumentException("THe matrix is not positive definite");
			}
			R[i][i] = sqrt(temp2);
			for(int j = i + 1; j < n; j++) {
				if (A[i][j] == 0.0) {
					R[i][j] = 0.0;
				} else {
					temp = 0.0;
					for(int k = 0; k <= i - 1; k++) {
						temp += R[k][i] * R[k][j];
					}
					R[i][j] = (A[i][j] - temp) / R[i][i];
				}
			}
		}
		return R;
	}
	
}